setwd("~/Dropbox/professional/Research/Active/Dyadic Design/")

set.seed(123)
library(sna)
library(ergm)

nsim <- 500

null.p <- numeric(nsim)
omit.p <- numeric(nsim)
miss.p <- numeric(nsim)

null.coef <- numeric(nsim)
omit.coef <- numeric(nsim)
miss.coef <- numeric(nsim)



for(i in 1:nsim){
	net <- rgraph(25,1,.1)
	net0 <- simulate(net~edges+istar(2)+istar(3), coef= c(-3.5,.75,-.1))
	x <- degree(net0,cmode="indegree")+ rnorm(25)*1.25
	dyadX <- t(net0[,])%*%net0[,]+matrix(rnorm(25*25,sd=3/4),25,25)
	null_est <- ergm(net0 ~ edges+edgecov(dyadX)+istar(2:3))
	miss_est <- ergm(net0 ~ edges+edgecov(dyadX)+gwesp(0,fixed=T)+gwodegree(1,fixed=T))
	omit_est <- ergm(net0 ~ edges+edgecov(dyadX))
	summary(omit_est)
	
	null.p[i] <- summary(null_est)$coefs[2,4]
	omit.p[i] <- summary(omit_est)$coefs[2,4]
	miss.p[i] <- summary(miss_est)$coefs[2,4]
	null.coef[i] <- summary(null_est)$coefs[2,1]
	omit.coef[i] <- summary(omit_est)$coefs[2,1]
	miss.coef[i] <- summary(miss_est)$coefs[2,1]
	print(i)
	
}

save(list=c("null.p","omit.p","miss.p","null.coef","omit.coef","miss.coef"),file="SimulationResults.RData")

# Type one error rate (two-tailed p-value)
p.vals <- seq(0.01,0.2,by=0.01)

null.error <- numeric(length(p.vals))
omit.error <- numeric(length(p.vals))
miss.error <- numeric(length(p.vals))
for(i in 1:length(p.vals)){
	null.error[i] <- mean(null.p < p.vals[i])
	omit.error[i] <- mean(omit.p < p.vals[i])
	miss.error[i] <- mean(miss.p < p.vals[i])
}

filei <- "type_one_error.pdf"
pdf(filei,height=3,width=3.5,pointsize=9,family="Times")
par(las=1,mar=c(4,5,1,1),cex.lab=1.25,cex.axis=1)
plot(p.vals,null.error,type="n",ylab="P(Type I Error)",xlab="",ylim=c(0,0.7))
abline(h=seq(.05,.7,by=.05),lty=3,col="grey70")
abline(v=seq(.025,2,by=.025),lty=3,col="grey70")
lines(p.vals,omit.error,lwd=1.5,col="red")
lines(p.vals,null.error,lwd=1.5,col="blue")
lines(p.vals,miss.error,lwd=1.5)
points(p.vals,omit.error,lwd=1.5,col="red",pch=4,cex=.85)
points(p.vals,null.error,lwd=1.5,col="blue",pch=1,cex=.85)
points(p.vals,miss.error,lwd=1.5,pch=6,cex=.85)
title(xlab="two-tailed p-value",line=2.25)
#abline(0,1)
legend("bottomright",legend=c("miss.pecified Model","Dyadic Independence","Full Model")[c(3,2,1)],col=c("black","red","blue")[c(3,2,1)],pch=c(6,4,1)[c(3,2,1)],bg="white")
dev.off()
